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ABSTRACT 

We consider the angular momentum exchange at the corotation resonance 
between a two-dimensional gaseous disk and a uniformly rotating external po- 
tential, assuming that the disk flow is adiabatic. We first consider the linear 
case for an isolated resonance, for which we give an expression of the corotation 
torque that involves the pressure perturbation, and which reduces to the usual 
dependence on the vortensity gradient in the limit of a cold disk. Although this 
expression requires the solution of the hydrodynamic equations, it provides some 
insight into the dynamics of the corotation region. In the general case, we find an 
additional dependence on the entropy gradient at corotation. This dependence 
is associated to the advection of entropy perturbations. These are not associ- 
ated to pressure perturbations. They remain confined to the corotation region, 
where they yield a singular contribution to the corotation torque. In a second 
part, we check our torque expression by means of customized two-dimensional 
hydrodynamical simulations. In a third part, we contemplate the case of a planet 
embedded in a Keplerian disk, assumed to be adiabatic. We find an excess of 
corotation torque that scales with the entropy gradient, and we check that the 
contribution of the entropy perturbation to the torque is in agreement with the 
expression obtained from the linear analysis. We finally discuss some implica- 
tions of the corotation torque expression for the migration of low mass planets 
in the regions of protoplanetary disks where the flow is radiatively inefficient on 
the timescale of the horseshoe U-turns. 

Subject headings: accretion, accretion disks — hydrodynamics — methods: nu- 
merical — planetary systems: formation — planetary systems: protoplanetary 
disks 
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1. Introduction 

It is known since the early eighties that low mass planetary objects (that is, up to a 
few Earth masses) embedded in protoplanetary gaseous disks should undergo a fast decay 
towards their central object, on timescales much shorter than the lifetime of the disk. This 
process, known as type I migration, has constituted for a long time a critical stage for 
the theory of giant planet formation. While it may account for the discovery of close-in 
extrasolar planets, with orbital periods of a few days, it renders problematic the build up 
of giant planet cores at distances of their central stars of several astronomical units. Most 
published studies of the tidal interaction of low mass objects with their parent disk have 
used either a barotropic assumption (such as a polytropic equation of state), or a locally 
isothermal equation of state. All these studies, whether analytical or numerical, confirmed 
the vigorous tidal interaction of the planet with the disk, leading to its inward migration on 
short timescales. 



There has been some exceptions to such assumptions: iMorohoshi &: Tanakal (120031 ) con- 
sidered the case of a planet interacting with an optically thin disk, in the shearing sheet 
approximation, and found that radiative effects can significantly alter the one-si ded torque 



between the planet and the disk. More recently, iPaardekooper fc Mellemal (120061 ) (hereafter 
PM06) have performed global, high resolution 3D calculations with nested grids that include 
radiative transfer. For the setup that they considered, they found that the total torque ex- 
erted by the disk on the planet increases with the disk opacity. For sufficiently large values of 
the opacity (and in the limit case of an adiabatic flow, corresponding to an infinite opacity), 
they find that the total torque on the planet is positive. This result is of great importance, as 
it potentially solves the lingering problem of type I migration. PM06 identified the existence 
of a hot, underdense part of the co-orbital region lagging the planet, which accounted for the 
torque excess that they measured. The present work corresponds to an attempt to further 
investigate this topic, so as to identify the physical mechanism responsible for these effects. 
For this purpose, we consider a more restricted situation, namely two-dimensional adiabatic 
fiows. 

This paper is organized as follows. In section [2] we set up the problem and define 
the notation. We then present an analysis of the corotation torque in an adiabatic disk 
in the linear regime, at an isolated resonance, at section [31 Our original motivation for 
the study of the linear regime was that PM06 found that the total torque reverses in a 
radiatively inefficient disk both for a 5 and a 0.5 planet, which pointed out that the 
effect is likely a linear one. In section HI we check by means of customized two-dimensional 
hydrodynamical simulations involving an isolated resonance the torque expression found in 
section [31 In section [SI we turn to the case of a planet embedded in an adiabatic disk, for 
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which we check that there is an excess of corotation torque that scales with the entropy 
gradient. We also check in this section that the torque excess corresponds to the sum of 
the linear contributions of all co-orbital corotation resonances, for a sufficiently small planet 
mass. We discuss the implications of the modified corotation torque expression for the issue 
of planet-disk tidal interactions, and we suggest further research on this topic in section [61 
We sum up our results in section [3 



2. Setup and notation 

We consider an inviscid, radiatively inefficient (that is to say, for our purposes, adiabatic) 
two-dimensional disk. In order to avoid corotation torque issues, we shall consider either 
a potential slowly turned on (sections [3] and H]) , or early stages after the introduction of a 
planet (section [5]) . The unperturbed state of the disk corresponds to a rotational equihbrium 
between the gravitational force of the central object, the pressure gradient and the centrifugal 
force. The unperturbed state is axisymmetric. The disk rotates with the angular speed fi(r), 
where r is the distance to the central object. We denote by p the pressure, E denotes the 
surface density, u and v respectively the radial and azimuthal velocities, ip the azimuthal 
angle. We denote by a "0" subscript the unperturbed quantities, and with a "1" subscript the 
perturbed ones. For instance, p(r, (yj) = poij-) +pi(r, y^). We shall essentially consider disks 
in which the unperturbed pressure and density are power laws of the radius, respectively 
with index A and cr: 

Po(^) oc r"'*' (1) 
So(r) cx r-\ (2) 

We shall make use of the two Oort's constants: 

. 1 rffi 

A = -r— , 3 
which scales with the local shear in the flow, and: 

B^L^-^^a + A, (4) 

2r dr 

which is half the vertical component of the flow vorticity, and which is also (2r)~^ times the 
radial derivative of the speciflc angular momentum. We will also use the epicyclic frequency 
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3. Linear analysis at an isolated resonance 

3.1. Basic equations 

We study the linear response of the disk to a perturbing non-axisymmetric potential 
$(r, yj) = cos[m{ip — Qpt)]. The perturbing potential rotates at constant angular 

velocity Qp. In the inertial frame, the linearized Euler equations of the disk are: 

^ + n^-2nv, = -^-^^ + ^^ (5) 

dt dip or So or Sg or 

and 

— -^ + 1]— -^H ui = — 7^ $ + — • (6) 

dt dip 2n rdip\ J:oJ 

The linearized continuity equation is: 

(9Si dEi 1 d , ^ , 1 d , 

dt dip r dr r dip 

We refer to the quantity S = pT,^'^ as the gas entropy, where 7 is the adiabatic index. 
The energy equation is equivalent in our case to the conservation of the gas entropy. The 
linearized conservation of the entropy along a fluid element path reads 

dSi dSi dSo 

where Si = S'o(pi/po~7Si/So)- We furthermore assume that the gas is described by an ideal 
equation of state so that po and Sq are connected by pq = SoC^/7, Cg being the adiabatic 
sound speed. 

We assume a perturbation of the form Xi^rn{f^) exp(im{(/? — fipt}) where Xi stands for any 
perturbed quantity of the fio\\ll]. We note Alu = m(Q,p — il) and we use the prime notation 
to denote d/dr. Eq. ([8]) can be recast as: 

E. = Pi + ^. (9) 

Combining Eqs. (jH]) and we are led to: 



— <^ $ + ^ ^} $ + ^ 

\l \ r \ r 



(10) 



We drop the subscript 771 ill 3^i^m(^) to improve legibility. 
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and 



K 



where S and V are given by 
and 

where \1/ is defined as 
and where JF is defined by 



s 



V 



1 dlnSo 
7 dlnr 

1 dlnpQ 
7 dlnr 

^ = pi/So 



with D = - Au"^ - SVcl/r"^. 

Substituting Eqs. (E]), ([10]) and dH]) into Eq. (ED leads to 

r2(<|. + m)" ^r{B + S){^ + - rS-^' + + = 0, 

where: 



i3= 1 + V 



dlnVt 
d\n.r ' 



C 



[r/S)' - 1 



V 



n 



Auj \ Auo^ 



and 



V 



(iln 
(ilnr 



(11) 
(12) 
(13) 
(14) 
(15) 



(16) 
(17) 



(18) 
(19) 
(20) 



Eq. (fT6|l reduces to the equation (15) of iLi et al.l (120001 ) if one considers the propagation 
of free waves ($ = 0), while it reduces to the equation (13) of iGoldreich fc Tremaind (119791 ) 
in the case of a homentropic {S = 0) fiow. 
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3.2. Corotation torque 



We now estimate the rate of angular momentum exchanged between the perturber and 
the radiatively inefficient disk described in section 13.11 This rate therefore corresponds to 
the disk torque, which we denote by F, and which we define as the torque exerted by the 
disk on the perturber (unless otherwise stated). It reads: 



f 9$ 
r = / Si(r, 09)— — r(ir(io9. 



(21) 



We limit ourselves to the torque exerted by the disk material lying in the vicinity of corota- 
tion, hence to the corotation torque, which we denote by Tc- In a linear analysis, this torque 
can be expressed as a series of contributions at each azimuthal wavenumber: Tc = J2m ^c,m- 
Each individual torque can be expressed, assuming that $ is real, as: 



(22) 



where denotes the imaginary part, is the corotation radius, and x = {r — rc)/rc. We 
substitute SqMi in Eq. ([9]) by the expression given bv Eg. (fTOj), and we k e ep on ly the terms 
which are large in the vicinity of corotation. As in lGoldreich fc Tremaind (119791 ). we assume 
that the disk responds to a slowly increasing perturbation and take Aa; to have a small, 
positive imaginary part a: 



where e 



Auj = m{Qp — Q) + ia ^ —mrfTt {r^[x + ie), (23) 
-a/[mrc^'{rc)\ > 0. In the vicinity of corotation, we can finally write: 



Ei(x) = ^(x) 



($ + ^)(x) 
X + ie 



2TS 



(24) 



We are primarily interested in the imaginary part of Si. In the limit e — > 0, we can write 
the terms that yield a non-vanishing contribution to the torque as: 



3 



[Si(a;)] = ^\^{x) 



-|- 7r5(x) 



X 



2TS 



(25) 



where S{x) is Dirac's delta function. The ffist two terms of the R.H.S. of Eq. (p5|) yield 
respectively the following contributions to the corotation torque: 



dx^[^{x)] 



2m7r2j^^5<l>($ + 3ft(^)) 



(26) 
(27) 
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The third term of Eq. fl25|) yields a contribution that can be shown to be neghgible, in 
the planetary context, compared to rc,m,2- This is shown in appendix Rl 

The first term, r^m,!; is the contribution of the function \E', such as in the barotropic 
case. The second term, Tc^rn,2, corresponds to a singularity at corotation, associated to a 
non- vanishing entropy gradient. It corresponds to the torque arising from the advection 
of entropy in the corotation region, which results in a surface density perturbation if the 
entropy is not uniform. The perturbation is singular for the surface density and the entropy, 
but not for the pressure (see section I3.2.2p . It remains confined to corotation, where it 
yields a singular contribution to the torque. Some further insight into the dynamics of this 
perturbation will be given in section 14. 3[ 

We provide in the next section an expression for the corotation torque in the limit of a 
cold disk, then we turn to the general case. 



3.2.1. Limit of a cold disk 

We contemplate here the case for which l^'l ^ |$|, which we shall refer to as a cold case. 
This condition depends on the strength of the perturbing potential, its radial scale, and on 
the disk temperature. In particular, in the planetary context, some corotation resonances 
may correspond to a cold situation, while others have |\E'| ~ |$|. Nevertheless, a given 
resonance eventually satisfies the cold case condition as the disk temperature tends to zero. 

The evaluation of Eq. fl26p requires an explicit expression for obtained by solv- 
ing the differential eq u ation (ITB]) in the vicinity of corotation. This has been done by 



Goldreich fc Tremaind (Il979l ) for a cold barotropic disk. An explicit solution can also be 
obtained for a cold adiabatic disk within the same level of approximation. Note however 
that some additional difficulties arise, in particular the existence of a double pole (term pro- 
portional to Auj^^) in the coefficients C and V, defined respectively by Eqs. ( ITSj) and (IT^ . 

We discard the double pole for the following reasons: 



• Unlike the simple pole, it scales with c^, which indicates that when the disk aspect 
ratio tends to zero, it becomes negligible; differently stated, there should be a critical 
disk thickness under which it is safe to neglect this term. 

• This term is the only one that depends both on the entropy and on the pressure 
gradients. As we shall see in section 15.2.21 our results of numerical simulations for a 
planet embedded in a disk with aspect ratio h = 0.05 show that the torque excess with 
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respect to an isothermal situation essentially depends on S, the gradient of entropy, 
which indicates that already for h = 0.05 the double pole term is negligible. 

The double pole is regularized with a very small amount of dissipation. Even the 
molecular viscosity suffices to render it negligible in the disks that we consider (S.-J. 
Paardekooper, private communication). 



Discar ding the double pole, and within the same level of approximation as lGoldreich fc Tremaine 
(Il979l ). Eq. (ITBl) can be recast, in the vicinity of the corotation, as 



where 



Pi 



(28) 



and q = {Dr/cs)r, ~ {i^r/cs). 



The general solution of Eq. (1281) reads 



Pi 



2q 



dt 



t + ie 



-qt 



dt 



t + ie 



(29) 



which reduces to the equation (53) of iGoldreich &: Tremaind (1l979l ) when S = 0. Combining 
Eqs. f l26|) and fl29|) yields the contribution Fcm,! to the corotation torque: 



Tc,m,l = To [{V + S) $2]^^ 



(30) 



where Tq = — (m7r^So)/(2i?ri7') is to be evaluated at the corotation radius. It can be 
approximated as (4m7r^So/3f2^)r^ in a Keplerian disk. 



The second contribution to the corotation torque, given by Eq. (1271) . is specific to the 
adiabatic case and involves the singularity arising from the entropy advection. Using Eq. ([2] 
and noting that |3fJ(\l')| ^ |$|, this contribution to the corotation torque reads 



c,ni,2 



(31) 



From Eqs. fl30!) and fl3T|) . we infer the corotation torque for a cold, adiabatic disk, which 
reads: 

r,,^ = ro[v<i>2]^^,^. (32) 

This expression does not depend on S. We note from Eqs. (|T5l) and fl20l) that V can be 
approximated as 

d\n So/5 



V 



dlnr 



(33) 
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since the disk aspect ratio at corotation h{rc) = c^(rr) /\rrfl(r r)] satisfies h(rr) <C 1. Eg. fl32l) 
therefore corresponds to the corotation torque expressioiu of Goldreich fc Tremaine (fl9^79i ). 
This argues that the corotation torque for a cold case does not depend on whether the 
disk can radiate energy efficiently (assuming a locally isothermal equation of state) or not 
(assuming an adiabatic energy equation). This can be expected on general grounds: in the 
cold disk limit, the internal energy of the fluid is negligible with respect to its mechanical 
energy, hence the power (and the torque) of the tid al force correspond to the cas e of non- 
interacting test particles, for which the expression of iGoldreich fc Tremaind (119791 ) prevails. 



3.2.2. General case 



We consider in this section the general case where we cannot neglect \E' with respect 
to $ in Eqs. (l26l) and (l27j) . as we have done in the previous section. Instead of re- 
sorting to a s olutio n of Eq. (ITB]) . we shall use a method similar to the method used by 



Tanaka et al. 



J2002h . based on the jump of angular momentum flux at corotation. In the 



case of Tanaka et al.l (I20021). t h is eve ntually yields a torque expression similar to the expres- 
sion of IGoldreich &: Tremaind (119791 ) . except that ^ has to be substituted by $ + (where 
7] is the enthalpy perturbation). The drawback of this method is that it provides a torque 
expression that depends on the (unknown) solution of the differential equation. Neverthe- 
less, it gives some insight into the dynamics of the corotation region, and allows to draw 
the general trends of the corotation torque in an adiabatic disk. In our case, the torque 
expression features \E' = pi/Sq. We note that in the isothermal case, IZhang fc Lail (120061 ) 
have provided an explicit solution for the perturbed enthalpy at corotation, that leads to a 
corotation torque expression that only depends on the forcing potential. 

We note that the jump of angular momentum flux at corotation misses the singular 
contribution of the entropy perturbation at corotation and as such leads only to an evaluation 
of Ec The contribution Ec,m,2 of the entropy perturbation to the corotation torque needs 
to be calculated similarly as in Eq. (IHTl) . The angular momentum flux Fa is given by: 



2n 

'0 



(34) 



where 3? stands for the real part and the star superscript denotes the complex conjugate. 
Eq. (I51|) can be written as Fa = J2m ^A,m with: 



FA,m = vrSor^ [^{uiMv^) + Q{u,)'^{v,)] . 



(35) 



^They have a negative sign because they consider the torque exerted by the perturber on the disk. 
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Combining Eqs. ffTOj) . ffTTl) and fl35l) . we obtain 

mvrEor 



(ir 



5 



-{3fJ($)^>(^) -9($)3fJ(^)} 
r 



(36) 



In the homentropic {S = 0) case, Eq. ( l36l) reduces to the expression used by iTanaka et al. 
( 12002| ). The contribution F^m,! to the corotation torque is then given by: 



c,m,l 



hm [Fa, 



(37) 



where > Vc and < rc are the radii of locations respectively beyond and before corota- 
tion, and where we evaluate the flux of advected angular momentum. 



Tanaka et al.l (120021 ) showed that $ + 77 is continuous at corotation. Here, since Eq. (fT6l) 



cannot be recast as an ordinary differential equation involving only ^-f^^, we have to consider 
more stringent albeit reasonable assumptions, namely that both $ and ^ are continuous at 
corotation. The fact that $ is continuous at corotation can be realized with an arbitrarily 
small softening length of the potential, in the case of an embedded point-like mass (for which 
the potential components would diverge logarithmically at corotation, in the absence of any 
softening). Assuming that $ is continuous at corotation, Eq. (1161) imposes that is also 
continuous at corotation (we would otherwise have a null linear combination of 6{x) and 
6'{x) functions with non- vanishing coefficients, which is impossible). 

Our continuity assumption implies that the terms proportional to S in the R.H.S. of 
Eq. (l36l) does not contribute to the torque. The jump in the advected flux therefore comes 
from the jump in (i($ + '^)/dr. 

We integrate Eq. (ITB]) over an infinitesimal interval containing r = r^. All finite terms 
in this equation yield a vanishing contribution, hence we are left only with the jump of 
(i($ + '^)/dr stemming from the second derivative term of Eq. f|T6|) and the poles of the 
terms C\E' and P$. This reads: 



d{^ + ^) 

I 

dr 



dr 



where 



2fi 



rVL 



7 (V + 2S) 



-[P.($ 



and Q 



vl/)(r,)-g$(r,)]. 



(3^ 



'2VL 



s 
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Using Eqs. ([36]), dST]), dMD and 25 = h."^ /2Vl, we find that 

Tc,m,i = To [{V + 25} 1$ + ^1' - 5 $ + ^ 
Eq. ( l39l) reduces to Eq. ( l30|) in the cold disk hmit. 



(39) 



We now come to the contribution Tc^m,2 of the entropy perturbation to the corotation 
torque. Eq. (!27|l yields: 

^e,„^,2 = -^o[5$3^^(<l> + ^)],^. (40) 

Eq. (HQ]) reduces to Eq. in the cold disk limit. 

The general expression for the corotation torque is obtained by accounting for the con- 
tribution given by Eq. fl39|) . and that of the entropy perturbation, given by Eq. fj40l) : 



To [{V + 25} 1$ + - 25 $ 3fJ(<l' + ^ 



(41) 



Eq. (HB reduces to the expression o f iTanaka et al.l (120021 ) when 5 = 0, while it reduces to 
that of iGoldreich &: Tremaind (119791 ) for a cold disk. 



A case of interest is that of a disk perturbed by a peaked potential (that of an embedded 
protoplanet for instance), for which |$+3?(\E^)| < |$|, and |<l>+3?(\l/)| < \'^{'^)\ at corotation. 
For such case, |$3?($ + \E')|r^ ^ |$ + ^E'l^^, hence the corotation torque may be approximated 
as Fc^m ~ — 2Fo[5$3?($ + '^)]rc- The corotation torque may therefore be much larger in 
the non-homentropic case (5 7^ 0) than in the homentropic case (5 = 0). Furthermore, its 
sign is given by that of 5 rather than that of V. This enhancement of the corotation torque 
in an adiabatic fiow may have a dramatic impact on the type I migration of an embedded 
protoplanet, as was noted by PM06. 



4. Numerical study of an isolated corotation resonance 

We check in this section the analytical predictions of section [3] by means of numerical 
simulations involving an isolated corotation resonance (hereafter CR). 



4.1. Numerical issues 



Our setup offers a number of similarities with the setup of iMasset fc Ogilvid (120041 ) 
for the case of an isothermal disk. The hydrodynamics equations for the disk described in 
section 13.11 are solved using the code Fargo. A description of the properties of this code 
is deferred to section [STT] in which the code is used to simulate an embedded planet. As in 
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Masset &: Ogilvid (120041 ). we deal with the m = 3 CR. The disk is therefore torqued by an 
m = 3 external potential $ that reads 



$(r, t) = T{t/T)<P{r) cos[3(^ - fi^t)], 



(42) 



where 0(r) denotes the radial profile of the potential, VLp its pattern speed (note that we 
work in the corotating frame), t is the time and where 



Tix) 



sin^(7rx/2) if a; < 1 
1 otherwise 



is a temporal tapering that turns on the potential on the timescale r. 

The total torque exerted by the disk on the perturber, given by Eq. fl?T]) . is evaluated 



by 



i=0 j=0 



2A(p 



(43) 



where Nr (Ng) is the radial (azimuthal) number of zones of the mesh, Sij is the surface area 
of zone ^ij and Sj ,,- are the external potential and surface density at the center of 

this zone, and A(f = 2tt/Ns is the azimuthal resolution. Furthermore, the contribution r^i 
of the function to the torque is obtained by substituting Ei by Pi/c^ in Eq. ([5T]). It is 
therefore evaluated by 



Tel 



i=0 j=0 



(44) 



where Pij and Csij are the pressure and sound speed at the center of zone (i, j). The 
contribution Tc,2 of the entropy perturbation to the torque is eventually estimated as follows: 



c,2 



Ec — Ec 1. 



(45) 



The radial computational domain is narrow enough to avoid the location of the m = 3 
inner and outer Lindblad resonances (see iMasset fc Ogilvid l2004l ). Despite this precaution, 
wave killin g zones next to the boun daries were implemented to minimize unphysical wave 
reflections ( de Val-Borro et al. 20061 ). Eurthermore, the torque evaluation is performed by 
summing only on a domain of the grid that does not contain the wave killing zones, and the 
summation includes a spatial tapering on the edges of that domain. 

The disk surface density and temperature are initially axisymmetric with power-law 
profiles: 

S(r) =E,(r/r,)-'^ (46) 
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and 

T{T)=T,{r/r:)-'+'f, (47) 

where and are the surface density and temperature at the corotation radius rc, and 
where / is the flaring index of the disk. The disk aspect ratio is given by h{r) = H{r)/r = 
h{rc){r /vc)^ , where H{r) is the disk scale height at radius r. A vanishing value of the flaring 
index / therefore corresponds to a uniform disk aspect ratio. The functions V and S are 
constant and read: 

V = 3/2 -a (48) 
S = a - (a + 1 - 2/)/7. (49) 



The main numerical parameters are those taken by lMasset fc Ogilvid ( 120041 ). namely a 
h{rc) = 0.01 disk aspect ratio at corotation, and Sc = 1. Our disk is inviscid. The libration 
islands are resolved by 30 zones azimuthally. As the potential increases, the radial width of 
the islands also increases. Their maximal radial width W spans approximately 20 zones. 

The results presented in next section have the following units: the mass of the central 
object M^, is the mass unit, the corotation radius Vc of our m = 3 CR is the distance unit 
and the Keplerian orbital period Torb at r = is 27i times the time unit. 



4.2. Results 

We consider three cases, corresponding respectively to Figs. [H [2^ and[2)D: 



1. An external potential with flat profile 0(r) = —10^^, as in iMasset fc Ogilvid (120041 ). 
This case, that we call the "flat potential case", has the following parameters: a = 2 
and / = —0.3, which implies, from Eqs. (jUj) and (H^ . that V = —0.5 and S ~ —0.57, 

2. A potential profile that corresponds to the m = 3 Fourier component of the smoothed 
potential of a M = 3.1 x 10~^M^, point-like object. The softening length is e = H{rc), 
which is approximately equal to W. The object rotates at speed Qp, with orbital 
radius r^. This neglects the pressure gradient effects, as we do not resolve the distance 
from orbit to corotation, but this distance is much smaller than the potential softening 
length, so this is not a concern in the present case. By contrast to the previous case, 
we call this situation the "peaked potential case" . The value of M was chosen so that 
0(^c) = — 10~^, as in the flat potential case. For this calculation we have a = 1.5 and 
/ = —0.3, so that V = and S ~ —0.71. The results are depicted in Fig. [2^. 
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ANALYTICAL 
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Fig. 1. — Corotation torque exerted by the disk on the perturber, as a function of time, 
assuming a flat radial profile of the potential. The results shown are obtained with an 
adiabatic calculation, except in the close-up, where we compare the isothermal and adia- 
batic corotation torques over the whole duration of the calculations. Numerical results are 
displayed with a symbol while the theoretical expectations are displayed with curves. We 
plot as a function of time the adiabatic corotation torque (diamonds and solid curve), the 
contribution of the function \E' to the torque (stars and dashed curve), and the contribu- 
tion of the entropy perturbation (triangles and dot-dashed curve). The long-dashed curve, 
which is nearly superimposed to the sohd curve, displays the corotation torque expression of 



Goldreich fc Tremaind (119791 ). The vertical solid line gives an estimate of the final libration 



time (see text). 
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3. A calculation similar to the previous one, except that a = 0.5 and / = —0.1, so that 
V = 1 and S ~ —0.71. The results are depicted in Fig. [2)d. 



For the three pairs (V,iS) quoted above, the tapering timescale value is r = 150Torb, 
which corresponds to the duration of the calculations. This is about three times larger than 
the final libration time, estimated as 

r,„~i(^)""^55r„,. (50) 

In each case we evaluate: 



• the total corotation torque (diamonds) with Eq. fH3|) . to be compared to the analytical 
expression (solid curve) given by Eq. fllT]) . In our units, Fq ~ 39.5, 

• the contribution of the function \E' to the corotation torque (stars) obtained with 
Eq. (jH]), the expected expression of which (dashed curve) is calculated using Eq. (!39|) . 

• the contribution of the entropy perturbation to the torque (triangles) using Eq. fH5|) . 
which is to compare to the prediction of Eq. fHOj) . represented by the dot-dashed curve. 



In these figures, the corotation torque first increases with time since the potential is 
progressively turned on until it reaches its final value at the end of the calculation. After 
some time it starts to oscillate. This oscillation corresp onds to the saturation of the CR, as 
the ratio t/Tub tends to unity ( jOgilvie fc Lubowll2003l ). Figs. [T] [2^ andEb therefore argue 
that our numerical simulations succeed in reproducing the results of our analytical study as 
long as t < Tiib, that is when a linear analysis is grounded (which requires that the time of 
the calculation be much smaller than the libration time). 

The examination of the results of these calculations leads to the following comments: 



In the flat potential case, depicted in Fig.dl we have 3fJ[V'(rc)] ~ — 0.02(/)(rc) throughout 
the calculation, where iplr) denotes the radial profile of This situation therefore 
corresponds to a cold case. As expected from Eg. (I32l) . the analytical corotation 
torque and the expression of iGoldreich fc Tremaind ( 19791 ) almost coincide. The close- 
up shows the torque evolution over the whole extent of the calculation, up to t = r. 
The torque obtained with a locally isothermal equation of state is also depicted. Our 
isothermal runs have same radial temperature dependence as the adiabatic runs (see 
Eq. p7|) . Although there is an entropy gradient in these isothermal calculations, it does 
not contribute to the corotation torque as it would in an adiabatic disk: the appearance 
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of the singular contribution at corotation in the adiabatic case is hnked (i) to the 
advection of entropy and (ii) to the appearance of a singularity in the perturbed density 
and temperature fields. In the isothermal situation, neither the entropy is conserved 
along a fluid element path, nor is a temperature singularity allowed to appear. The 
comparison of isothermal and adiabatic calculations shows that, as expected for a cold 
case, the adiabatic and isothermal torques coincide, as long as we are in the linear 
regime. We note that both torques do not oscillate about since the potential reaches 
a stationary value only at the end of the calculation. 

For the two calculations of the peaked potential case, depicted in Figs. [2^ and [2)d, we 
find that 'Si[ip{rc)] ~ -O.20(rc). Thus, the term | - $3?($ + slightly dominates the 
term |$ + \E'p in Eq. fl41l) . Because S < for these calculations, the corotation torque 
in the adiabatic case (diamonds and solid curve) is larger than the corotation torque 
in an isothermal dis k (long dashed curve) with the same parameters, as predicted by 



Tanaka et al.l (|2002| ). In particular, in the case for which V = 0, the isothermal coro- 
tation torque vanishes, while we find a net, positive corotation torque for an adiabatic 
flow, in correct agreement with the analytical expression. 



4.3. Dynamics of the corotation region 

We discuss in this section the dynamics of the corotation resonance of an adiabatic disk 
and give some comments about the corotation torque expression of Eq. fHTl) . 

In the isothermal case, the corotation torque expression involves the p roduct of the 



gradie nt of vortensity and the square of the effective potential ($ + 77), see e.g. iTanaka et al. 



( I2OO2I ). The torque is then given by the angular momentum budget between material flowing 
outwards and material flowing inwards at corotation, regardless of the sign of $ + 77. Eq. 0411) 
displays a term that has a similar behavior, except that it does not feature the vortensity 
gradient only, but rather V + 2S. This factor sca les with the (l ogarithm ic) gradient of 



(En/-B )5'^/'^, which is a key quantity considered by iLi et al.l (I2OOOI ) and by iLovelace et al. 



( 119991 ). who pointed out that vortensity is not conserved in a two-dimensional adiabatic flow. 



In addition to this term, Eq. fHT]) contains a term that scales with $[$ + 9ft(\I')]. The 
sign of this term therefore depends on the relative signs of $ and $ + 3?(\I'). In order to get 
some insight into the physical meaning of this term, we show at Fig. [3] the response of the 
disk in the corotation region, for the entropy and the surface density. These flelds correspond 
to the calculation with the flat potential proflle considered at the previous section. The disk 
has a negative radial entropy gradient. Therefore, libration brings the (larger) inner entropy 
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Fig. 2. — Same as Fig. [H except that the results are obtained with a peaked potential, 
with V = ( V = 1) in the left (right) panel. The long-d ashed curve in both panels shows 
the expectation from the corotation torque expression of iTanaka et al.l (120021 ). denoted by 
TTW02. 




Fig. 3. — Relative perturbation of entropy (left) and surface density (right) for an isolated 
resonance, at t ~ 1.5 Tub. Libration is clockwise. 
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to the outer part of the hbration islands, yielding a positive perturbed entropy (brighter 
areas), while it brings the (smaller) outer entropy to the inner part of the libration islands, 
yielding a negative perturbed entropy (darker areas). An opposite behavior is observed for 
the perturbed density, since the relative pressure perturbation (not represented) is much 
smaller. 

The sign of this torque component can be understood as follows. Fig. H] depicts the 
situation in two cases: $ and $ + K(\E') have same sign (left), and $ and $ + 3?(\Ef) have 
opposite signs (right). In the left case, the negative perturbed surface density on the outside 
of corotation is located in the region where dcp^ < 0, hence the perturbation yields a positive 
torque on the perturber. A similar conclusion applies to the material flowing inwards which 
has positive perturbation of surface density. The torque on the perturber is therefore positive, 
in agreement with the sign of — + 3fJ(\E')]. An opposite conclusion holds for the case 
where $[$ + ^{^)] < 0. 

The order of magnitude and functional dependence of this torque component can be 
justified as follows. As the sign has been justified at the previous paragraph, we give here an 
estimate of the absolute value. The perturbed surface density on the outside of corotation 
is ~ |iSSo5/rc|, where 5 = [($ + 3fJ(\E'))/(— 8^41?)]^/^ is an order of magnitude of the width 
of the libration islands. The specific torque in the region of surface density perturbation 
is ~ |m$|, while the area covered by the perturbation of surface density scales with 
The torque arising from this region therefore scales with |mrc5^iS$So|, which is exactly the 
scaling of IFq + 5R(\1')]|, within a numerical factor in 0{1). 

The singular behavior of this torque component, which stems from Eq. (HOi) . and which 
appears as a Dirac's delta function at corotation, can be understood as follows: as the 
strength of the perturbation decreases, the width of the libration islands tends to zero, 
while the libration time tends to infinity (libration disappears), hence we are left, in the 
linear regime, with a torque contribution that comes strictly from the corotation radius and 
therefore appears as singular. 

It is worth noting that only half of the second term of Eq. fHTl) comes from Eq. fHU]) . 
Eq. fl5Ul) . which is obtained from the momentum flux jump, and which as such captures effects 
occurring at a finite (albeit small) distance from corotation, also displays a term similar to 
that of Eq. (l40l) . The advect ion of entropy perturbati ons is not a silent process: it triggers the 
emission of pressure waves ( Foglizzo fc Taggeilboool ). Our torque expression indicates that 



half of the energy required to advect entropy in the libration islands is evacuated through 
pressure waves. 
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5. Application to the case of an embedded protoplanet 

In section |3], we derived an expression for the corotation torque between a radiatively 
inefficient disk and an external rotating potential. This expression is successfully reproduced 
by local numerical simulations of an isolated corotation resonance, in the linear regime. 
We now contemplate the case of an embedded protoplanet in a radiatively inefficient two- 
dimensional disk, for which all co-orbital corotation resonances are simultaneously active. 



5.1. Numerical features and setup 



Our numerical simulations are performed with the code Fargo. It is a staggered mesh 
hydrocode that solves the Navier-Stokes, continuity and energy equations on a p olar grid. 
It us es an upwind transport scheme with a harmonic, second-order slope limiter (Ivan Leer 
19771 ). Its particularity is to use a change of rotating fram e on each ring of the polar grid, 
which increases the timestep sig nificantly JMassetl EoOOal jbl) . thereby lowering the compu- 
tational cost of a given calculation. The energy equation that we implemented in Fargo 



is: 



de 
dt 



+ V.(ev) = -pV.v + g, 



(51) 



where e is the thermal energy density, v = {u, rQ)^ denotes the flow velocity, p is the 
vertically integrat ed pressure and Q is a heating source term that accounts for the disk 
viscosity (see e.g . D'Aneelo et al.ll2003l ). The energy equation solver is implemented as in 



Stone fc NormanI fll992[ ) 



In this work, the disk is taken inviscid so Q = 0. There is no radiative transfer either, 
since the disk is assumed to be radiatively inefficient. Furthermore, p and e are connected 
by an ideal equation of state p = (7 — l)e, where the adiabatic index 7 is set to 1.4. This 
equation of state can be expressed in terms of the disk temperature T and surface density 
E as p = ST. The adiabatic sound speed reads Cs = VlT, hence Cg = ^/yCg^iso, where c^^iso 
refers t o the isothe rmal sound speed. We comment that the Lindblad torque, which scales 
(jWardlll997l ). is therefore weakened by a factor of 7 in an adiabatic disk. The same 



as c„ 



is true of the corotation torque, when there is no entropy gradient. We checked both effects 
with appropriate calculations, not reproduced here. This plays in favor of a total torque 
reversal in adiabatic disks with a negative entropy gradient. 

The disk is initially slightly sub-Keplerian (the pressure gradient is accounted for in 
the centrifugal balance), axisymmetric, with power-law profiles for the surface density and 
temperature given by Eqs. (H6l) and (H71) . 
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For a comparative purpose, calculations involving a locally isothermal equation of state 
are performed. In isothermal calculations, no energy equation is solved: p and E are simply 
connected hj p = Sc^ j^^. These isothermal calculations have same initial surface density and 
temperature profiles as the adiabatic runs. 

The disk is perturbed by the smoothed potential of a protoplanet. We adopt a Plummer 
potential, with a softening length e = 0.6H{rp) (unless otherwise stated), being the planet 
orbital radius. This fiducial value is quite substantial for our purposes, but investigating the 
disk response at much smaller softening lengths, where the adiabatic effects on the corota- 
tion torque are increasingly important, requires a very large resolution. A high resolution 
systematic study at small softening length will be presented in a forthcoming work. 

The protoplanet is held on a fixed circular orbit, at r = r^. The disk parameters are 
summed up in Table [H where they are expressed in the following unit system: is the 
length unit, the mass of the central object is the mass unit and (GM*/rp^)~^/^ is the 
time unit, G being the gravitational constant (G = 1 in our unit system). We denote by Torb 
the planet orbital period, Mp the planet mass and q = Mp/M^ the planet to primary mass 
ratio. 



5.2. Results 

5.2.1. An illustrative example 

We show the results of an illustrative 
primary mass ratio (corresponding to Mp = 
The horseshoe libration time is 

Tlib = 

Table 1. Reference parameters. The disk is inviscid 



Parameter 


Notation 


Reference value 


Aspect ratio at r = Vp . . 


Kvp) 


0.05 


Surface density at r = rp 




2 X 10-^ 


Softening length 


£ 


0.03 


Adiabatic index 


7 


1.4 


Mesh inner radius 


''^min 


0.4 


Mesh outer radius 




1.8 


Radial zones number .... 


Nr 


512 


Azimuthal zones number 


Ns 


2048 



calculation with ag = 2.2xl0 planet to 
7.3 if the central object has a solar mass). 
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where is the protoplanet angular velocity and Xs denotes the half-width of the horseshoe 



region. iMasset et al.l (120061 ) have given an estimate of Xg in the isothermal case, that reads 
Xg ~ 1.16rp^^yq/h{rp). A streamline analysis was performed and confirmed that this estimate 
holds for an adiabatic disk, if one substitutes h{rp) with ^h{rp). We find therefore rub ~ 
60 Torb- Numerical diffusion eventually alters the conservation of entropy. Nevertheless, the 
horseshoe region spans 20 zones radially, which is sufficient to follow the horseshoe dynamics 
over several libration times. Since we are concerned here with a fraction of the libration 
time, the entropy is conserved with a good level of accuracy over the duration of our runs, 
and it can be regarded as a Lagrangian tracer of the flow. 

Two calculations were performed: an adiabatic and an isothermal one. Both lasted 
thirty orbital periods, hence half the horseshoe libration time. This calculation has a = 0.5 
and / = 0, as in PM06. This gives S ^ -0.57. 

Fig. [5] displays the gas entropy, surface density and pressure obtained in the adiabatic 
calculation, after 15 Torb- Each field represents the relative perturbation of the corresponding 
quantity with respect to the unperturbed state. For instance, the top right panel shows 
[S(r, — So(r)]/So(r). While the azimuthal range spans the whole [0, 27r] interval, the radial 
range depicted is restricted to a band of width 2.5xs around the corotation radius r^. We 
overplot streamlines to the entropy panel to give an idea of the extent of the horseshoe region. 
The vertical dashed line represents the corotation radius. Whereas the pressure panel does 
not display any significant perturbation, the entropy and density panels show the propagation 
of a perturbation inside the horseshoe region, which slides along the separatrices. This is 
reminiscent of the behavior commented in the case of an isolated resonance at section 14.31 

The interpretation of this dynamics is as follows: the entropy of the fiuid elements is 
conserved as they perform a horseshoe U-turn in the co-orbital region. When there is initially 
an entropy gradient at corotation, the co-orbital dynamics yields an entropy perturbation 
that has a sign opposite of that of the entropy gradient on the outwards U-turns, and the 
sign of the entropy gradient on the inwards U-turns. Since the pressure field is only weakly 
perturbed, the entropy perturbation is related to a density perturbation of opposite sign 
and, in relative value, of same order of magnitude. Therefore, if there is a negative entropy 
gradient at corotation {S < 0, as in the example shown here), the co-orbital dynamics yields 
a negative density perturbation at ip < ipp and a positive density perturbation at ip > ipp, 
with straightforward consequences for the corotation torque. Using an expression inherited 
from the terminology of Riemann solvers, we call this perturbation a contact discontinuity. 
A contact discontinuity is characterized by a discontinuity in the density and temperature 
fields, while the pressure and velocity fields are continuous. A contact discontinuity is simply 
advected by the fiow. Here it follows the horseshoe dynamics, and it remains confined to the 
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horseshoe region. 

We give hereafter a simple estimate of the relative perturbation of the disk surface den- 
sity due to the advection of entropy. We consider a fluid element that performs a horseshoe 
U-turn from the inner part of the horseshoe region (where we assume that there is no entropy 
perturbation, which is true as long as t < Tiib/2) to the outer part. All physical quantities 
at the inner (outer) leg of the horseshoe streamline are denoted by a minus (plus) subscript. 
A first-order expansion yields, assuming no pressure perturbation: 

P± = Poircjil =F Ax/rJ, (53) 

where < x < is the distance of the streamline to corotation, and: 

S_ = So(r,)(l + ax/r,). (54) 

On the outer horseshoe leg, the disk surface density is perturbed according to the entropy 
perturbation and reads: 

S+ = So(r,)(l + i?-ax/r,), (55) 

where R is the relative perturbation of surface density aX r = rc + x (we assume a symmetric 
horseshoe U-turn), due to the entropy advection. Entropy conservation along the fluid 
element path {S^ = S^) leads to: 

R = 2-(a--] =2-S. (56) 
Tc \ 7/ 

The horseshoe U-turn that we have considered lags the planet ((/? < ipp). A similar conclusion 
holds for a horseshoe U-turn that switches from the outer leg to the inner one (at ip > (pp), 
hence we finally have: 

R{x) = 2xS/rc, Vx G [-x^, +Xs]. (57) 

The bottom right panel of Fig. [5] displays the slices of the perturbed density field at t = 
15Torb, for (f — ipp = 1 (diamonds) and ip — ipp = —1 (stars). The two horizontal dashed lines 
display the values of R{—Xs) and R{xs), where x^ is estimated through a streamline analysis. 
Similarly, the long-dashed curve shows R{x) = 2xS/rc, which is in correct agreement with 
the calculation results. The surface density structure in the horseshoe region is therefore 
dictated by the sign of S. In particular, we do not expect any contact discontinuity in the 
homentropic case (5 = 0). We have checked this prediction with a numerical simulation (not 
presented here). 
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5.2.2. Excess of corotation torque and entropy gradient 
An order of magnitude of the excess of corotation torque arising from the perturbation 



(Ward 


1991; 


Masset 


2001) 



i?(x)So(rc). We consider the outwards horseshoe U-turns that occur ai ip < ipp. Assuming, 
in this order of magnitude estimate, that the rotation profile of the disk is unperturbed, we 
evaluate the variation of angular momentum flux of the horseshoe disk material after the 
U-turn attributable to the change of the disk's surface density: 

ArHS-= I \-2Ax)^^R{x){3, + 2Br,x)dx, (58) 
Jo 

where jc is the specific angular momentum of the material at corotation. The first factor of 
the integrand of Eq. flSSl) represents the material velocity in the corotating frame, due to the 
shear. The last factor is the material specific angular momentum obtained from a first order 
expansion at corotation. Similarly, we obtain the change of angular momentum flux due to 
the perturbation of surface density on inwards horseshoe U-turns: 



ATus+ = {-2Ax)^QR{-x){ic-2Br^x)dx. (59) 
Jo 

Adding Eqs. (158!) and (!59|) . we are left with: 

AFhs = 2 / {-2Ax) ■ T.oR{x) ■ 2Brcxdx = -AABEoSxt- (60) 
Jo 

Fig. E] shows the excess of corotation torque between an adiabatic and isothermal calculation 
with same parameters, as a function of the half-width of the horseshoe region. This excess is 
obtained by subtracting the total torque of an adiabatic and an isothermal calculation (the 
isothermal torque being rescaled by a factor 7"^, since Cs = Y^c^.iso)- We call this difference 
the torque excess for further reference. Each data point corresponds to a calculation with a 
given planet mass, for which we determine Xg through a streamline analysis. We find that 
the torque excess approximately scales as x^, and that it is within a factor 2 of our order of 
magnitude estimate, given by — AFhs- 

The torque expression of Eq. (HTl) as well as the horseshoe drag expression of Eq. (!60l) 
suggest that the torque excess scales with S, hence with the entropy gradient. In order to 
check that, we have undertaken a number of calculations with different values of S. These 
calculations have q = 2.2 x 10~^, and the disk parameters are those of Table [TJ Each entropy 
gradient is realized with different combinations of the indexes of the pressure and surface 
density power laws. Adiabatic effects on the torque are assessed in two different ways: 
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1. By calculating the torque excess, as in Fig. [61 



2. By evaluating the following integral: 




clJ d<^ 




(61) 



which provides an estimate of the torque due to the contact discontinuity (this contri- 
bution arises from perturbations of S which do not have a pressure counterpart). In 
the linear regime, Eq. fl6Tl) amounts to a summation over m of the last term of Eq. (124]) . 
We shall check this statement in the next section. 

These two estimates of adiabatic effects on the torque value are shown respectively in 
Figs. [7K and [7)d. Remarkably, they coincide within ~ 25 %. We will comment further this 
coincidence in the next section. 

The main conclusion that can be drawn from the results of Fig. [7] is that the torque excess 
(or the contact discontinuity contribution) essentially depends on the entropy gradient, as 
expected. The excess is positive for a negative entropy gradient, hence we may expect the 
total torque exerted on a planet embedded in a radiatively inefficient disk to be a positive 
quantity if the radial entropy gradient is sufficiently negative. 



We have given at Eq. fHOj) an estimate of the singular torque contribution from the 
contact discontinuity at an isolated resonance, while we have estimated the total contribution 
in the planetary case of the contact discontinuity using Eq. fl^ at section 15.2.21 We check 
in the present section that this total contribution corresponds to the sum over m of the 
torque expression of Eq. (HOj) . For this purpose, we have adopted a planet to primary mass 
ratio g = 5 X 10~^, as the one adopted in the previous sections {q = 2.2 x 10^^) led to 
poor agreement, presumably because of the onset of non-linear effects. For each azimuthal 
wavenumber m, we measure 3?(\Efm) from the calculation output (at t = 5Torb), and we 
evaluate the sum over m of the torque Tc^rn,2- 



5.2.3. Connection to the analytical expression 



r, 



oo 



lim T' 



(62) 



fc— >+oo 



where: 




m<k 



(63) 
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is the partial sum of Tc,rn,2- We compare the torque contribution given by Eq. fl6T]) to Too- 
The results are presented in Fig. [HI The agreement between the direct torque measurement 
and the linear estimate is excellent. 

This confirms what we anticipated in section [3l2l and what is shown in appendix [XJ that 
the contribution of the last term of Eq. (!25|) to the torque is negligible in the planetary con- 
text. Also of interest is the torque density associated respectively to p/c^ and S —p/cj.. The 
sum of these two torque densities is the total torque density. They are represented at Fig. O 
The total torque density displays a smooth profile and a narrow peak at corotation. This is 



remin iscent of the torque density found by PM06 (their Fig. 2) or by iMorohoshi fc Tanaka 



(I2OO3I ) (their Fig. 3). The decomposition above sphts this total torque density in a smooth 
component arising from p/c^, which reminds the torque density in an isothermal disk, and 
a sharp, localized torque density arising from S — p/c^. This corresponds to the torque 
density of the contact discontinuity contribution given by Eq. (1611) . Fig. [9] shows that this 
contribution (which is singular at corotation in the linear case for an isolated resonance) is 
here bounded by the extent of the horseshoe region. 

We comment the surprising agreement found at the previous section between the torque 
excess and the contribution of the contact discontinuity. The linear analysis suggests that 
the former should be the sum of Fq,™ [25|$m + \E'mp - 25$m($m + 3fJ(\&m))], which, in 
the limit where |$m + 3fJ(^m)| < \^{^m)\ and |$m + 3fJ(^m)| < should reduce 

to —2TQmS^rn{^m + ^{"^ m)) , that is twice the contribution of the contact discontinuity 
(see section 13.2.21) . Nevertheless, for the calculations presented here, we can check that 
2 J2m ''^l^m + is almost exactly compensated by m^rni^m + 5^(^m))- Namely the 
ratio of the former to the latter quantity is found to be 1.07, which explains why the full 
excess expression essentially amounts to the contact discontinuity contribution. Presumably 
this coincidence is fortuitous and linked to the relatively large softening length that we use. 
At smaller softening length, the term in $m(^m + 3^(^m)) should largely dominate, yielding 
a ratio of 2 between the torque excess and the contribution of the contact discontinuity. We 
note that PM06 also quote that the torque estimate given by their equation (1) accounts for 
the total torque within 25 % (this equation can also be seen as an estimate of the contact 
discontinuity contribution). This seems to suggest that the softening length of 0.6H{rp) 
that we adopted is a correct choice to reproduce the magnitude of the corotational effects in 
adiabatic three-dimensional disks. 
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6. Discussion 
6.1. Softening length 



In an isothermal disk, the corotation torque scales with |$ + r^p (ITanaka et al.ll2002l ). 
Even if $ diverges at corotation, $ + remains finite, which makes the isothermal corotation 
torque a quantity relatively insensitive to the softening length. The situation is quite different 
for the effects linked to the entropy advection that we present here: they involve the product 
$($ + \[^), which diverges when $ does. Adiabatic effects on the corotation torque should 
acquire a very large magnitude at small softening length. In particular, if the softening length 
is smaller than the distance from orbit to corotation, the magnitude of these effects should 
strongly depend on this distance, which scales with the pressure gradient. If one regards 
the softening length as a proxy for the altitude in a three-dimensional disk, the extent of 
the disk vertical scale height concerned by these very small softening length issues should be 
small, however, since the distance from orbit to corotation is a fraction of Tph?. Nevertheless, 
it is of interest to investigate the behavior of the corotation torque in an adiabatic fiow at 
very small softening length to assess the importance of such effects. Owing to the very large 
resolution required to investigate this problem, we defer this investigation to a forthcoming 
work. 



6.2. Saturation 



The origin of the effects presented here is the advection of entropy in the corotation 
region, that triggers an entropy perturbation (and therefore a density perturbation) when- 
ever there is an entropy gradient in the equilibrium profile. Libration occurs on different 
timescales for the different streamlines of the corotation region, which tends to stir the en- 
tropy and to flatten out the entropy profile across the corotation region (be it the horseshoe 
region in the planetary case or a libration island in the isolated resonance case). This is 
quite similar to the behavior of the corotation torque in an isothermal disk, which tends 
to saturate because the vortensity profile is fiattened out by libration. In this case, it is 
the viscous diffusion which can prevent the flattening out of the proflle if it acts sufliciently 
rapidly to establish the large scale gradients b e fore a libr ation time. This h a s been studied 
for an isolated resonance bv IColdreich fc Saril J2003h and lOgilvie fc Lubowl (120031 ) and by 
Balmforth fc KorycanskyI (120011 ) and iMassetl (120011 ) for a planetary co-orbital region. In 
both cases, the degree of saturation of the corotation torque in steady state depends on the 
ratio of the libration time and of the viscous time across the libration region. The dissi- 
pative processes required to prevent the torque saturation in the situation presented here 



-27- 



should be able to impose the large scale entropy gradient over the corotation region in less 
than a libration time. Radiative processes (cooling and heating) should therefore occur on 
a timescale longer than a horseshoe U-turn (otherwise the flow can rather be considered 
as locally isothermal), but they should act on a timescale shorter than the libration time. 
We provide an estimate of the horseshoe U-turn time and of the libration time for a small 
mass object embedded in a gaseous disk. The horseshoe half-width Xs is ~ ^p^J q/h{rp). 
Neglecting pressure effects and writing a simplified Jacobi constant for a test particle near a 
horseshoe U-turn as: J = —GMp/ {2Brp\ip — ipp\) + A{r — r^)^, we can estimate the distance 
of closest approach between the planet and a test particle flowing along a horseshoe separa- 
trix as rp|Ay9|s = VL^H{rp) /2\AB\ = 0{H{rp)). The time required to perform a horseshoe 
U-turn can be deduced using the radial drift velocity of the test particle when it crosses the 
orbit, at its closest approach from the planet: x = GMp/ {2BrpA<f1) . That yields: 

ru^turn = 2xs/x = Q^p h{rpf' q''/'/ {A' B) 

where Rh = rp{q/3y^^ is the Hill radius of the planet, and where the last equality holds 
for a Keplerian disk. When the planet emerges from the disk {H{rp) ~ Rh), the horseshoe 
U-turn occurs on the dynamical timescale. When dealing with an embedded object however 
{Rh < H{rp)), the horseshoe U-turn time can be substantially longer than the dynamical 
time (e.g. 10 times longer for an Earth mass object embedded in a disk with h{rp) = 0.05). 

Using Eq. (l52l) . we are led to: 

^ h{rp)-\ (65) 

''"U— turn 

There is at least an order of magnitude difference between the horseshoe U-turn time and the 
libration time in a thin disk, hence it should be possible to find a location in the disk where 
the cooling time is much longer than the U-turn time and yet shorter than the libration 
time, so as to maintain an unsaturated corotation torque. 
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6.3. Entropy gradient and baroclinic instability 

The effect that we present in this two-dimensional analysis is of particular interest 
when there is a negative entropy gradient at corotation, since this may suffice to halt 
type I migration. It would be of interest to generalize the present analysis to the case 
of a three-dimensional baroclinic disk. We comment also that in such syste ms, a nega- 



tive entropy gradient may render the disk unstable to a baroclinic instability (IKlahrl 12004 
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Klahr &: Bodenheimeiil2003l ). It is certainly important to examine the interplay of the baro- 



clinic instability and of the corotational effects presented here. The turbulence generated by 
the baroclinic instability, in particular, could provide a mechanism to prevent the saturation 
of the corotation torque, much like the turbulence arising from the MRI can prevent the 
corotation torque saturation in an isothermal disk. 



7. Conclusions 



We evaluate the corotation torque between an adiabatic gaseous disk and a uniformly 
rotating external potential. In the linear case for an isolated resonance, we find a singular 
contribution at corotation which scales with the entropy gradient, and which arises from the 
advection of entropy within the libration region. This effect neither exists in isothermal or 
locally isothermal flows, nor does it exist for barotropic fluids (such as fluids described by a 
polytropic equation of state). We provide a torque expression at an isolated resonance which 
involves the pressure perturbation at corotation. We then check the torque expression by 
two-dimensional adiabatic calculations that involve an isolated resonance. In particular, we 
exhibit a case with a fiat vortensity profile, for which the corotation torque does not cancel 
out and is in correct agreement with the analytical expression. We then turn to the case of 
an embedded planet, for which we find an excess of corotation torque in the adiabatic case, 
which scales with the entropy gradient. For a sufficiently small planet mass, we check that 
this excess can be accounted for by a summation over the resonances of the torque excess 
that we found in the first part. This confirms that this effect is essentially a linear effect. 
We finally discuss in section [6] some open questions linked to the softening length, to the 
saturation, to the case of a three-dimensional baroclinic disk, and to the interplay with the 
baroclinic instability, on to which theoretical efforts should focus in a nearby future. 



We wish to thank Sijme-Jan Paardekooper and John CB. Papaloizou for interesting 
discussions on the topics covered in this manuscript. We also thank Alessandro Morbidelli 
for a thorough reading of a first version of this manuscript, and an anonymous referee for 
comments that led to an improvement of the paper. 



A. Additional contribution to the corotation torque 

We check hereafter that the contribution of the last term of Eq. fl2S]) is negligible for an 
embedded planet. For this purpose, we compare G = J dx'^[^{x)]^{x)/x to — 7r[3fJ(\Ef)$],.^. 
Fig. [TOk shows the m = 8 component of \E' = p/So for the calculation presented in sec- 
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tion I5.2.3[ We clearly see that the behavior of 53(\E') in the vicinity of corotation comes 
from the overlap of the behavior arising at the inner and outer Lindblad resonances. Be- 
tween its outermost inner minimum at r_ ~ 0.87 and its innermost outer maximum at 
r+ ~ 1.12, 53(\E') can be considered as having a linear dependence in x. Notwithstanding 
the decrease of $ as one recedes from corotation, the main contribution to G will come 
from Q'[\Ef(x)]/x between these two radii, as it exhibits a flat behavior over this range. This 
yields G ~ 2$(rc)|^j[\E'(r-|-)]|, that it to say a result comparable in order of magnitude to 
— 7r[3?(\E')$]j.^. Nevertheless, the final contribution of the extra term is much smaller than 
the singular one at corotation for the following reasons: 



1. The potential decreases sharply as one recedes from corotation, which provides a cut-off 
to the extra term, that is not localized at corotation. 

2. The extra term is partially compensated for by the first term of Eq. ffTOl) . which we 
have neglected in writing Eq. fl2^ . and which yields another term in Eq. fl25|) that 
reads —{J-'S/r'^Q)r^d'^[^{x)]/dx. Adding this additional term and the last term of 
Eq. (^B^, we are left with: 



ci3[^(a;) 
dx 
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which shows that, in addition to the potential cutoff, the last term of Eq. (l25l) should 
be decreased by a factor of 4 (we neglect, at this level of accuracy, the jump in 53(\1'') 
at corotation). 



We have checked on the calculation presented in section 15.2.31 that the contribution of these 
extra terms is indeed small compared to the singular contribution at corotation at all m. 
This is shown in Fig. [TOb . from which we can conclude that the total contribution of the 
extra terms is about an order of magnitude smaller than the singular contribution. The 
agreement that we found in section 15.2.31 between the numerical simulation of an embedded 
planet and the torque series, which was of the order of a percent, might then be fortuitous. 
Nevertheless, we expect an agreement of the order of 10 %, still very satisfactory. These 
findings are also compatible with the fact that we hardly see any diffuse torque density 
outside of the horseshoe region in Fig. [9]d. 
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Fig. 4. — Sketch of the corotation region when $ and $ + 3fJ(^) are in phase (left) and when 
$ and $ + 3?(^) are in opposition (right). The minima and maxima of $ are indicated at the 
left, while the minima and maxima of $ + 3?(^) are indicated at the right. In the corotation 
region, material hbrates about the maxima of the effective potential $ + 3?(^'). We assume a 

negative entropy gradient, hence material flowing outwards has a negative perturbed surface 
density, while material flowing inwards has a positive perturbed surface density, as indicated 
by the minus and plus signs. 
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Fig. 5. — Top row and bottom left: relative perturbations of the gas entropy, surface density 
and pressure, at t = 15 Torb ~ Tiib/4. The protoplanet is located in r = Vp, ip = ipp. In the top 
left panel, streamlines are overplotted and the vertical dashed line stands for the corotation 
radius Tc- In the top right and bottom left panels, the color scale is adjusted to highlight the 
advection of the entropy perturbation (see text). The nearly horizontal overdensity structure 
ai if = ifp is the protoplanet 's wake. Bottom right: slices of the relative perturbed density 
field at the same time, a.t (p — ipp = 1 (diamonds) and ^ — ^p = —1 (stars). The two horizontal 
dashed lines refer to the values of R{—Xs) and R{xs), while the long-dashed curve displays 
the quantity 2(r — rc)S /tc (see text and Eq. (!57j) ). 
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Fig. 6. — Torque excess (see text) as a function of the half-width of the horseshoe region. 



2 o.in-n> 



o 





^ 0=0.5 X=1. 5 










A a=^.0\=2^ 










W o=1.5X=2.9 










0=0.5 1=1 
«=1-5>^== 


.1 
.5 










^ a=O.OX=0.0 










X (i=1.0it=1.4 










O o=1.5X=2.1 












V <i=0.5X=0.3 










> <j=1.5X=1.7 












0=1.0 X=0.7 








S 


0=0.5 i=0.0 










o=1.5 X=1.4 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 

S = CT - Wy 



4-1 0" 



-2.10" 



o=1.5X=2.9 



0=0.5 X=1.1 
c=1.5X=2.5 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 

S = CT - Wy 



Fig. 7. — Torque excess (left) and contact discontinuity contribution to the torque (right) 
as a function of S. Although the calculations display some scatter for a given value of «S, 
the different points can be considered as aligned within a good level of approximation. The 
slope of the dependence is negative. 
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Fig. 8. — Partial sums of the torque series given by Eq. (jHSD (diamonds) and direct estimate 
of the contact discontinuity contribution, given by Eq. (!6T|) (dashed hne). The asymptotic 
value of the partial sum almost coincides with the direct estimate (i.e. the diamonds almost 
lie on the dashed line at large m), hence with a very good accuracy we have Too = ^cd (see 
text). 



- 36 - 




1.04 1.06 



Fig. 9. — Left: total torque density (solid curve) and torque density of p/c^ (dashed curve). 
Right: torque density of E — p/c^. The vertical dashed line shows the corotation radius, 
while the two vertical dotted lines show the extent of the horseshoe region. 
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Fig. 10. — Left: m = 8 azimuthal component of = p/Sq (the real part is shown by a 
solid line, the imaginary part is shown by a dashed line). The dotted line shows the m = 8 
component of the potential (which is purely real). Right: singular contribution of at 
corotation (solid line), contribution of the extra term in 55[\E'(x)]/x of Eq. (^31) (dotted line), 
contribution of the first term of Eq. ffTOl) (dashed line), and total contribution of these extra 
terms (dash-dotted line). Each contribution is normalized to the maximum value of the 
singular contribution of ^ at corotation. 



